/* 
 * File:   main.cpp
 * Author: paweln66
 *
 * Created on 27 październik 2012, 17:35
 */

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <time.h>
#include <cmath>

/*
 *grain representation
 */
class Grain{
public:
    int grainID;
    double ro;
    bool recrystalized;
    bool recInPrvStep;  //TO REMOVE

    Grain(){
        grainID = 0;
        ro = 0;
        recrystalized = false;
        recInPrvStep = false;
    }
};

/*
 *This class is used to run this model
 */
class ModelRunner {
private:
    Grain ** space;
    Grain ** spacePrv;

    int dimX, dimY;
    int numGrains;
    int progress;   //progress to finish naiv grain growth
    int numSteps;   //number of steps in process
    double ro_Intial, ro_0, ro_critical;
    double A, B;
    double t,  delta_t;

    /*Variable for GPGPU computing*/
    int * cl_spaceGrainID;  //cl_uint
    int * cl_spaceGrainIDPrv;   //cl_uint
    float * cl_spaceGrainRo;    //cl_float
    float * cl_spaceGrainRoPrv;    //cl_float
    int * cl_spaceGrainRecryst;     //cl_uint
    int * cl_spaceGrainRecrystPrv;     //cl_uint
    int cl_width, cl_height;  //cl_uint
    float cl_deltaRo, cl_Rocritical;   //cl_float

    /*Private methods*/
    int periodicMinusX(int index);
    int periodicMinusY(int index);
    int periodicPlusX(int index);
    int periodicPlusY(int index);
    int naiveGrainGrowth();
    int recrystalizeGrain(int i, int j);
    void setGrainMoore(int i, int j);   //TODO templatre for this method which can be use for all neighbourhoods
    void setGrainNeumann(int i, int j);
    void updateSpace();
    void drawInTerminal();  //TO remove
    void drawInFile(int s = 0);  //TO remove
    void setDeployment();
    void updateDeplyment(int i, int j, double t);
    void startRDX();
    bool inicializeSpace();
    bool isOnBoardNeumann(int i, int j);
    double calculateDeployment(double t);

public:
    /*COnstructors and destructors*/
    ModelRunner(){
        space = 0;
        numGrains = 0;
        progress = 0;
    }
    ModelRunner(int DimX, int DimY, int NoGrains){
        this->numGrains = NoGrains;
        this->progress = 0;
        this->dimX = DimX;
        this->dimY = DimY;
        this->t = 0;
        this->delta_t = 0.001;
        this->A = 86710969050178.5;
        this->B = 9.41268203527779;
        this->ro_critical = 4215840142323.42;

        space = new Grain * [dimX];
        spacePrv = new Grain * [dimX];
        for (int i =0; i< dimX; i++){
            space[i] = new Grain [dimY];
            spacePrv[i] = new Grain [dimY];
        }
    }

    ~ModelRunner(){
        if(space != 0){
            for(int i =0; i<dimX; i++)  delete [] space[i];
            delete [] space;
        }
        if(spacePrv != 0){
            for(int i =0; i<dimX; i++)  delete [] spacePrv[i];
            delete [] spacePrv;
        }
    }

    /*Public methods*/
    bool runModel(int NoSteps);
};

/*
 * 
 */
int main(int argc, char** argv) {

    ModelRunner runner(10, 10, 4);

    //std::cout << "Check: " << runner.start() << std::endl;
    runner.runModel(10000);
    
    return 0;
}

/*Method to run all recrystallization process*/
bool ModelRunner::runModel(int NoSteps){
    this->t = 0;
    this->numSteps = NoSteps;

    this->inicializeSpace();
    this->drawInFile();
    this->drawInTerminal();
    this->naiveGrainGrowth();
    this->drawInTerminal();
    this->startRDX();
    this->drawInTerminal();

    return true;
}

int ModelRunner::naiveGrainGrowth(){
    //for( int step = this->progress; step <= this->dimX * this->dimY; step++){
    while ( this->progress < this->dimX * this->dimY){
        for (int i = 0; i<dimX; i++){
            for( int j =0; j<dimY; j++){
                this->setGrainNeumann(i, j);
            }
            this->drawInFile();
        }
        updateSpace();
        this->drawInFile();
    }

    return 0;
}

bool ModelRunner::inicializeSpace(){
    srand ( time(NULL) );

    int pX, pY;

    /*for(int i = 1; i <= this->numGrains; i++){
        pX = rand() % this->dimX;
        pY = rand() % this->dimY;

        space[pX][pY].grainID = i;  //only one problem will be when both random number will be the same
        spacePrv[pX][pY].grainID = i;

        progress++;
    }*/

    //TO REMOVE
    progress = this->numGrains;
    space[1][1].grainID = 1;
    spacePrv[1][1].grainID = 1;
    space[3][1].grainID = 2;
    spacePrv[3][1].grainID = 2;
    space[4][3].grainID = 3;
    spacePrv[4][3].grainID = 3;
    space[5][2].grainID = 4;
    spacePrv[5][2].grainID = 4;

    return true;
}

void ModelRunner::setGrainNeumann(int i, int j){
    if (space[i][j].grainID != 0) return;   //if grain is setted do nth

    int nghID[4];
    /*Vertical*/
    nghID[0] = spacePrv[i][periodicMinusY(j)].grainID; //up
    nghID[1] = spacePrv[i][periodicPlusY(j)].grainID;  //down
    /*Horizontal*/
    nghID[2] = spacePrv[periodicMinusX(i)][j].grainID; //left
    nghID[3] = spacePrv[periodicPlusX(i)][j].grainID;  //right

    //set grain for max number of neighbour
    int * grainNghCount = new int [this->numGrains];
    for (int g = 0; g<numGrains; g++) grainNghCount[g] = 0;
    //count neighbour for all grains
    int notNgh = 0;
    for (int cell = 0; cell<4; cell++){
        if(nghID[cell] == 0)    notNgh++;
        for (int grain = 1; grain<=numGrains; grain++)
            if (nghID[cell] == grain) grainNghCount[grain - 1]++;
    }
    if (notNgh == 4)    return;     //if grain doesn't have neighbour
    //find id grains with max No of neighbours
    int * maxNghID = new int [this->numGrains + 1];   //it should be a list, but right now I dont know if I can use SDT library
    maxNghID[0] = 0;  //ID grains with max number of neighbours, grain id with max number ngh will on the top of the list
    int max = grainNghCount[0];    //amount of neighbour
    int numMaxNgh = 1; //how many is maximum neighbours grain ID
    for (int g = 1; g <= numGrains; g++) {
        if(grainNghCount[g] > max){
            numMaxNgh = 1;
            max = grainNghCount[g];
            maxNghID[0] = g;
        }
        else if(grainNghCount[g] == max){  //if more than one neighbour
            numMaxNgh ++;
            maxNghID[numMaxNgh - 1] = g;
        }
    }

    //set id grain(i, j) with max number of neighbour
    if (numMaxNgh == 1) space[i][j].grainID = maxNghID[0] + 1;
    else {
        space[i][j].grainID = maxNghID[rand() % numMaxNgh] + 1;
    }
    this->progress++;

    delete [] maxNghID;
    delete [] grainNghCount;
}

int ModelRunner::periodicMinusX(int index){
    index = index -1;
    if(index < 0) return this->dimX - 1;
    else return index;
}

int ModelRunner::periodicMinusY(int index){
    index = index -1;
    if(index < 0) return this->dimY - 1;
    else return index;
}

int ModelRunner::periodicPlusX(int index){
    index = index + 1;
    if(index > this->dimX - 1) return 0;
    else return index;
}

int ModelRunner::periodicPlusY(int index){
    index = index + 1;
    if(index > this->dimY - 1) return 0;
    else return index;
}

void ModelRunner::drawInTerminal(){
    for (int i = 0; i < dimX; i++){
        for (int j = 0; j < dimY; j++)
            std::cout << space[i][j].grainID << " ";
        std::cout << std::endl;
    }

    std::cout << std::endl;
}

void ModelRunner::drawInFile(int s){
    std::ofstream file;

    file.open("output.txt", std::ios::app);

    file << "Step: " << s << std::endl;
    for (int i = 0; i < dimX; i++){
        for (int j = 0; j < dimY; j++)
            file << space[i][j].grainID << " ";
        file << std::endl;
    }

    file << std::endl;

    file.close();
}

/*
 *TO copy spaca to spacePrv
 */
void ModelRunner::updateSpace(){
    for (int i = 0; i < dimX; i++)
        for (int j =0; j<dimY; j++){
            spacePrv[i][j].grainID = space[i][j].grainID;
            spacePrv[i][j].recrystalized = space[i][j].recrystalized;
            spacePrv[i][j].ro = space[i][j].ro;
            spacePrv[i][j].recInPrvStep = space[i][j].recInPrvStep;
        }
}

double ModelRunner::calculateDeployment(double t){
    double d1, d2, d3;
    d1 = A/B;
    d2 = 1 - A/B;
    d3 = std::exp((-B) * t);
    return d2 * d3  + d1;
}

/*
 *set initial deployment
 */
void ModelRunner::setDeployment(){
    this->t += this->delta_t;
    this->ro_Intial = calculateDeployment(t);
    this->ro_0 = this->ro_Intial;

    for (int i = 0; i<this->dimX; i++)
        for (int j = 0; j<this->dimY; j++){
            space[i][j].ro = ro_Intial;
            spacePrv[i][j].ro = ro_Intial;
        }
}

/*
 *dyslocation gun to set assign deployment (only for next space)
 */
void ModelRunner::updateDeplyment(int i, int j, double deltaRo){
    

//    double ro1 = calculateDeployment(t);
//    double deltaRo = ro1 - this->ro_0;
//    this->ro_0 = ro1;

    int rand = std::rand() % 4;

    if (this->isOnBoardNeumann(i, j)){
        double delta = (1.0 - (double)rand * 0.1) * deltaRo;
        space[i][j].ro += delta;
        spacePrv[i][j].ro += delta;
    }
    else{
        double delta = (double)rand * 0.1 * deltaRo;
        space[i][j].ro += delta;
        spacePrv[i][j].ro += delta;
    }
}

/*
 *verify if grain is on board
 */
bool ModelRunner::isOnBoardNeumann(int i, int j){
    int activeGrainID = spacePrv[i][j].grainID;

    /*Vertical*/
    if (spacePrv[i][periodicMinusY(j)].grainID != activeGrainID && !spacePrv[i][periodicMinusY(j)].recrystalized) return true; //up
    if (spacePrv[i][periodicPlusY(j)].grainID != activeGrainID && !spacePrv[i][periodicPlusY(j)].recrystalized) return true;  //down
    /*Horizontal*/
    if (spacePrv[periodicMinusX(i)][j].grainID != activeGrainID && !spacePrv[periodicMinusX(i)][j].recrystalized) return true; //left
    if (spacePrv[periodicPlusX(i)][j].grainID != activeGrainID && !spacePrv[periodicPlusX(i)][j].recrystalized) return true;  //right

    return false;
}

/*
 * verify if neighmbour grain is recrystalized and return grain id if yes or 0 if not
 */
int ModelRunner::recrystalizeGrain(int i, int j){

    int nghID[4];
    int numRecrystGrain = 0;    //number of recrystalized grain neighbour
    /*Vertical*/
    if(spacePrv[i][periodicMinusY(j)].recrystalized){
        nghID[numRecrystGrain] = spacePrv[i][periodicMinusY(j)].grainID;  //up
        numRecrystGrain++;
    }
    else if(spacePrv[i][periodicPlusY(j)].recrystalized) {
        nghID[numRecrystGrain] = spacePrv[i][periodicPlusY(j)].grainID;  //down
        numRecrystGrain++;
    }
    /*Horizontal*/
    else if(spacePrv[periodicMinusX(i)][j].recrystalized) {
        nghID[numRecrystGrain] = spacePrv[periodicMinusX(i)][j].grainID;  //left
        numRecrystGrain++;
    }
    else if(spacePrv[periodicPlusX(i)][j].recrystalized) {
        nghID[numRecrystGrain] = spacePrv[periodicPlusX(i)][j].grainID;   //right
        numRecrystGrain++;
    }
    else {
        return 0;
    }

    if (numRecrystGrain == 1) return nghID[0];
    else return nghID[rand() % numRecrystGrain];
}

/*
 *run dynamic recrystallization process
 */
void ModelRunner::startRDX(){
    this->setDeployment();
    this->progress = 0;

    for (int step = 0; step < this->numSteps && this->progress < this->dimX * this->dimY ; step ++){
        this->t += this->delta_t;
        double ro_1 = calculateDeployment(t);
        double deltaRo = ro_1 - this->ro_0;
        this->ro_0 = ro_1;

        for (int i = 0; i < this->dimX; i ++){
            for (int j = 0; j < this->dimY; j ++){
                this->updateDeplyment(i, j, deltaRo);

                space[i][j].recInPrvStep = false;

                /*IF grain is on board and ro > ro_critical*/
                /*new embryo*/
                if (spacePrv[i][j].ro > this->ro_critical && this->isOnBoardNeumann(i, j) && !spacePrv[i][j].recrystalized){
                    space[i][j].recrystalized = true;
                    space[i][j].grainID = ++this->numGrains; //new embryo grain
                    space[i][j].ro = ro_Intial;
                    space[i][j].recInPrvStep = true;
                    this->progress++;
                }

                /*IF neighbor was recrystalized in previous step (t-1) and grain[i, j].ro > ro_in_all_ngh*/
                else if (/*spacePrv[i][j].recInPrvStep*/ !spacePrv[i][j].recrystalized &&
                        spacePrv[i][periodicMinusY(j)].ro < spacePrv[i][j].ro &&    //up
                        spacePrv[i][periodicPlusY(j)].ro < spacePrv[i][j].ro &&   //down
                        spacePrv[periodicMinusX(i)][j].ro < spacePrv[i][j].ro &&  //left
                        spacePrv[periodicPlusX(i)][j].ro < spacePrv[i][j].ro)  //right
                {
                    int newNgh = this->recrystalizeGrain(i, j);
                    if(newNgh != 0){ 
                        space[i][j].grainID = newNgh;
                        space[i][j].recrystalized = true;
                        space[i][j].ro = ro_Intial;
                        space[i][j].recInPrvStep = true;
                        this->progress++;
                    }
                }
            }
        }
        this->updateSpace();
        this->drawInFile(step);
    }
}